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Abstract 

We present an update of the CLUMPY code for the calculation of the astrophysical /-factors (from dark matter annihilation/decay) for 
any Galactic or extragalactic dark matter halo including substructures: halo-to-halo concentration scatter may now be enabled, boost 
factors can include several levels of substructures, and triaxiality is a new option for dark matter haloes. This new version takes 
advantage of the cf itsio and HEALPix libraries to propose fits output maps using the HEALPix pixelisation scheme. Skymaps 
for y-ray and v signals from generic annihilation/decay spectra are now direct outputs of CLUMPY. Making use of HEALPix routines, 
smoothing by a user-defined instrumental Gaussian beam and computing the angular power spectrum of the maps is now possible. 
In addition to these improvements, the main novelty is the implementation of a Jeans analysis module, to obtain dark matter density 
profiles from kinematic data in relaxed spherical systems (e.g., dwarf spheroidal galaxies). The code is also interfaced with the 
GreAT toolkit designed for Markov Chain Monte Carlo analyses, from which probability density functions and credible intervals 
can be obtained for velocity dispersions, dark matter profiles, and /-factors. 
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PROGRAM SUMMARY 

Program Title: CLUMPY 
Programming language: C/C+-I- 
Computer: PC and Mac 
Operating system: UNIX(Linux), MacOS X 

RAM: between 500MB and 1GB depending on the size of the 
requested skymap 

Keywords: dark matter, indirect detection, Jeans analysis, y-rays, v 
Classification: 1.1, 1.7, 1.9 

External routines/libraries: CERN ROOT (http: //root. corn. ch), 
GSL (http://www.gnu.org/software/gsl), cfitsio (http: 
//heasarc.gsfc.nasa.gov/fitsio/fitsio.html), HEALPix 
C-l—I- and F90 (http://healpix.sourceforge.net/index.php), 
GreAT (http://lpsc.in2p3.fr/great) (for MCMC analyses 
only), and Doxygen (http: //www. doxygen. org) (optional) 

Nature of problem: Calculation of dark matter profile from kinematic 
data, y-ray and v signals from dark matter annihilation/decay. 

Solution method: Solve the integro-differential Jeans equation 
(optimised for speed) for several generic distributions (dark matter 
profile, light profile, velocity anisotropy). Integration of the DM 
density (squared) along a line of sight for generic dark matter haloes 
with substructures (spatial, mass, concentration distributions). Draw 
full skymaps of y-ray and v emission from dark matter structures, 
smoothed by an instrument PSF using HEALPix tools. 

Restrictions: The diffuse extragalactic contribution to the signal (and 
y-ray attenuation) as well as secondary radiation from dark matter 


remain to be included in order to provide a comprehensive description 
of the expected signal. 

Running time: This is highly dependent of the user-defined choices of 
DM profiles, precision e and integration angle crint: 

• ~ 1 hour for a full skymap (including substructures) with ofint = 
O.r and6 = 0.01; 

• < 1 mn for a 5° x5° skymap (including substructures) with = 
O.r and6 = 0.01; 

• ~ 5 mn for a typical Jeans/MCMC analysis (on a ‘ultrafaint’-like 
dwarf spheroidal galaxy) using a constant anisotropy profile. 
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1. Introduction 

Dark matter indirect detection aims at measuring the end 
products of dark matter (DM) annihilation or decay p, d, y, 
y). The recent results from the PAMELA [1,2] and AMS-02 [3] 
experiments for charged particles, and the wealth of Fermi-LAT 
results for y-rays [4, 5,6], show that we are starting to probe the 
region of interest of the parameter space for new physics. Due 
to the complexity of the signal and backgrounds involved, the 
need for public tools for cross-checks and progresses in the field 
is mandatory. 

Several particle physics public tools exist to calculate the 
spectrum of species produced from dark matter annihilation and 
decay (e.g., micrOMEGAs [7, 8], DarkSUSY [9]). The astrophys- 
ical side of the calculation depends on the nature of the created 
particle: for charged species p, d), a diffusion/convection 
equation must be solved, and several propagation codes are 
available (e.g., GalProp^ and Dragon^). Neutral particles (y, 
y) propagate in straight lines, and the main uncertainty in the 
signal is related to the /-factor calculation from DM structures 
and substructures: CLUMPY^ [10] (Paper I) is the only public 
code for the full /-factor calculation. 

CLUMPY has been used to calculate /-factor in DM haloes of 
dwarf spheroidal galaxies (dSphs) [11, 12], DM-supported Hi 
clouds [13], and towards the Galactic centre to place limits on 
DM annihilation using the ANTARES neutrino telescope [14]. 
The code has been extended to provide additional outputs (sort¬ 
ing, population study) when applied on a sizeable sample of 
galaxy cluster haloes [15, 16, 17], or galactic subhaloes. The 
second CLUMPY release, presented in this paper, includes these 
additional outputs, but also provides a better description of sev¬ 
eral quantities related to DM haloes and their substructures, as 
well as more fiexible and useful outputs. The main novelties of 
the code are the following: 

1. implementation of the well-established Jeans analysis [18] 
to extract DM profiles from the kinematics of baryonic 
tracers. It is interfaced with a Markov Chain Monte Carlo 
(MCMC) engine which enables the user to perform the 
complete analysis from kinematic data to the astrophysi- 
cal factors probability density distributions; 

2. triaxiality of the DM haloes (Milky Way, or any other DM 
halo) can now be enabled; 

3. use of the HEALPix (Hierarchical Equal Area isoLatitude 
Pixelation) and cfitsio input/output libraries for 2D 
skymaps fits outputs, including smoothing by a Gaus¬ 
sian beam and power spectrum tools (the code for 2D 
skymaps has also been optimised for speed); 

4. use of tabulated DM y-ray and y spectra from [19]^ to 
compute fiuxes for annihilation and decay (y oscillations 
are included). 


^http://galprop.Stanford.edu 

^http://www.dragonproj ect.org 

^http://Ipse.in2p3.fr/clumpy 

^http://www.marcocirelli.net/PPPC4DMID.html 


This paper highlights CLUMPY’s main features, old and new; 
a more thorough description (and examples) can be found in 
the fully Doxygen-documented code. The paper is organised as 
follows. Section 2 briefiy recalls the formalism for calculating 
the /-factor in the presence of dark matter substructures. Sec¬ 
tion 3 provides updated formulae for this new version, which 
includes a statistical description of the concentration param¬ 
eter, the multi-level boost calculation, and triaxiality for DM 
haloes. Section 4 focuses on the description of the new outputs 
provided in this version (fits files containing maps using the 
HEALPix pixelisation scheme, y-ray and y fiuxes). Section 5 
deals with the Jeans analysis, its implementation in CLUMPY, 
and the interface with the MCMC engine GreAT. Run exam¬ 
ples with a short description of CLUMPY (and its installation) are 
given in Section 6. We conclude in Section 7. 


2. Dark matter annihilation or decay: reminder 

2.7. Fluxes 

The y-ray or y fiux dOy^y/dE" from dark matter annihilat¬ 
ing/decaying particles is expressed as the product of a particle 
physics term by an astrophysical contribution /. At energy E 
and in the direction {if/, 0), the fiux integrated over the solid an¬ 
gle AQ = 2;: (1 - cos Ofint) is given by 

d<Dyv 

-^(E, 4,, e, AQ) = X J{4,, e, AQ), (1) 

aE aE 

in which dQ = dyS sin ada is the elementary solid angle around 
the line-of-sight direction ij/, 6 (longitude and latitude in Galac¬ 
tic coordinates)^. 


2.7.7. Particle physics term 

The particle physics term depends on whether the DM can¬ 
didate annihilates or decays. In this version (as in the previous 
one), we only consider the continuum emission [e.g., 19]: 


(^ann^) 


do. 


dE 




dN: 


y,y 


dE 


Bf X 




TDM ^DM 


(annihilation) 


(decay) 


( 2 ) 

with moM the mass of the DM candidate, Bf the branching ratio 
into the final state / and its yield per reaction dNy^yjdE (see 
§4.1), and 


• CTann IS the annihilation cross section, and {cr^^^v) the an¬ 
nihilation rate averaged over the DM velocity distribution, 
6 equals 2 (resp. 4) for a Majorana (resp. Dirac) fermion; 

• tdm is the decay lifetime. 


^ A sketch of the coordinate framework is provided in the code documenta¬ 
tion (see also Fig. 6 of Paper I). 


2 











2.1.2. Astrophysical J (annihilation) or D (decay) factor 

The astrophysical factor relies on the integration over the 
solid angle AQ of some power of the DM density p(</r, 6, /, a,{3) 
at coordinate (/, a^jl) in the line-of-sight direction 0)\ 

/ r ir 1 ^ f (annihilation) 

d/dQx ^ ^ (3) 

Jo Ji.o.s [p (decay) . 

Note that, depending on the community, J- and D-factors 
are either expressed in astrophysics units (Mq kpc“^ and 
Mq kpc“^ respectively) or particle physics units (GeV^ cm“^ 
and GeV cm“^ respectively). All calculations in CLUMPY are 
performed in astrophysics units, but a new keyword intro¬ 
duced in this version (gSIMU_IS_ASTRO_OR_PP_UNITS) allows 
the user to select the preferred units for the outputs (plots, 
ASCII and fits files). 

In the following, we concentrate on the /-factor calculation 
for annihilation, for which the contribution from substructures 
is able to boost the signal (there is no boost for decaying DM). 
In both cases, the formalism is the same and is implemented 
similarly in CLUMPY. 

2.2. J-factor and substructures: formalism 

Here, we briefly recap the formalism used to take into ac¬ 
count substructures in CLUMPY. We refer the reader to Paper I 
for a more detailed description. The changes brought by this 
new release regarding the handling of substructures are post¬ 
poned to §3. 


The calculation of /subs and /cross-prod described in Paper I 
was using only one level of substructures. This new release of 
the code allows the inclusion of more levels of substructures as 
described in §3.1. 

2.2.2. Continuum limit 

A typical DM halo of lO^^M© (Milky-Way like) contains up 
to 10^"^ substructures, which renders the explicit calculation of 
the above sums prohibitive. This huge number allows the use 
of the continuum limit as the clump positions and masses are 
random variables, drawn from distribution functions obtained 
by N-body numerical simulations and/or semi-analytical calcu¬ 
lations^. 

As detailed in §3.1, the above equations can often be replaced 
by averaged ones. In particular, the total averaged density cor¬ 
responds to 

(Ptot)(^) — Psm(^) + (Psubs)(^)9 (8) 

where Psm(^) is the smooth component and (psub)(^) the spher¬ 
ical shell average density of substructures (at each radius r). 
The total averaged spherical shell density <ptot)(^) is usually 
parametrised by Zhao or Einasto profiles (see §3.3.4 of Paper I 
and prof iles . h). A saturation density psat (gDM_RH0SAT) pro¬ 
vides a cut-off radius below which the annihilation rate is con¬ 
stant (equilibrium between free fall time and annihilation time). 
For instance, selecting the total average density profile (ptot) 
and the clump distribution parameters (psubs) (see next para¬ 
graph), Eq. (8) is used to get the smooth distribution Psm(^) that 
is plugged in Eq. (5). 


2.2.1. Formal description 

From Eq. (3), the total DM density p must be known at 
each position. In the ACDM cosmological model, structures 
form in a bottom-up manner: micro-haloes form first, larger 
ones collapse later, and this process, accompanied with merger 
events, lead to the global picture of clumps within clumps 
within clumps, etc. Each DM clump can be seen as a density 
peak inside its host halo, and it is therefore convenient to sepa¬ 
rate, for a given halo, the main distribution (called smooth halo) 
from the contributions of each clump. 

The astrophysical contribution to the annihilation flux is thus 
explicitly written to be 


J = 



where p^j corresponds to the inner density of the i-th clump con¬ 
tained in the volume element. Three terms arise from this equa¬ 
tion (smooth only, clumps contribution, and cross-product): 


J subs — 


/c 


cross-prod 



plmsLX 

PsmYpi^ldQ. (7) 

0 //min 


2.2.3. Validity of the mean description 

At first order, a random variable (e.g. the mass and position 
of substructures) is described by its average value and variance. 
Departure from the average can arise if a small number of ob¬ 
jects contribute significantly to the total /-factor, which hap¬ 
pens if a massive clump dominates, or if one of the smallest 
halo (the smaller, the more numerous they are) is sitting almost 
at the observer location. The latter configuration only happens 
for clumps in the Galaxy, since clumps in dSphs or extragalactic 
objects are far away. 

As presented in Paper I, for /-factor skymaps of the Galaxy, 
CLUMPY relies on a combination of the calculation of the aver¬ 
age signal and the calculation of individual drawn clumps above 
and below a critical distance /cut (the computation of which is 
detailed in Paper I^), respectively. This strategy ensures a con¬ 
trolled and extremely quick calculation of skymaps: the number 
of clumps to draw in the Galaxy is reduced from a few tens of 
thousands to a few hundreds (see table 1 of Paper I) depend¬ 
ing on the user-chosen level of precision (or more precisely, the 
level of fiuctuation selected w.r.t. the mean signal). For all other 


^Large uncertainties remain on these distributions, all the more because 
small halo masses are not resolved, even in the most computationally heavy 
simulations: CLUMPY is partly designed to enable quick calculations of the J- 
factor for any input distribution, in order to check the sensitivity/robustness of 
the results against the uncertain parameters of the distributions. 

^The critical distance is obtained by requiring the relative error of the signal 
integrated from /cut to remain lower than a user-defined precision requirement. 
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objects—like dSphs, galaxies, galaxy clusters—, the mean de¬ 
scription is usually sufficient. 

3. Updates and novelties in the /-factor computation 

Three major improvements have been included: the concen¬ 
tration of the clumps is now dealt with a distribution function, 
so as to include an uncertainty around the mass-concentration 
relationship, and the calculation of the boost due to substruc¬ 
tures can include, if required, several levels of substructures 
(§3.1); it is now possible to deal with triaxial DM haloes in 
addition to spherical ones (§3.2). 

3.1. Improved concentration description and subhalo levels 

For a given DM profile, the physical properties of a subhalo 
are fully defined by its position, mass, and concentration^. In 
the previous CLUMPY release, the position and mass were ran¬ 
dom variables of user-defined distribution functions. The mass- 
concentration relation was fixed so that two DM haloes of equal 
mass would have had the same concentration. Numerical simu¬ 
lations have however shown significant uncertainties in the de¬ 
termination of the mass-concentration relation, which could be 
parametrised by a dispersion around an average relation. There¬ 
fore, in this release, the concentration parameter is a new ran¬ 
dom variable, characterised by a specific distribution function 
as described below. 


The latter case means that the concentration of a clump 
of mass M is randomly drawn from the above dis¬ 
tribution around the mean concentration c(M) (calcu¬ 
lated from a gENUM_CVIRMVIR enumerator), with a scat¬ 
ter (Tc{M) (e.g., 0.14 - 0.18 [20, 21, 22]). The 
parameters in clumpy_params.txt (see table 2) are 
gDM_FLAG_CVIR_DIST and gDM_LOGCVIR_STDDEV. 

3.1.2. CLUKPY formulae to calculate (/) 

The description of the concentration in terms of a distribu¬ 
tion function implies some modifications to the average (/) 
as defined in paper I (see §3.3.6 and the documentation in 
clumps. h for more details). The formulae below recap and 
extend Eqs. (17) - (22) of Paper I. The average mass and the 
mean of (some power of) the distance are left unchanged and 
read 


<M): 


= / M^dM, <f>=| ( r^Cdidn. 

-'Mmin dM Jo J/^i„ dV 


( 11 ) 


The intrinsic luminosity of a clump J1(M, c) now depends on 
the concentration c. 


£(M,c)=f pl(M,c)dV, (12) 

-'Vcl 

SO that the mean luminosity (over a range of M and c) becomes 


3.1.1. Substructures: random variables and distributions 
For a total number of clumps A/tot in a host halo, the substruc¬ 
ture distribution is modelled by: 


d^N 

dVdMdc 


dPv, . . 

dV dM dc 


(9) 


In the above equation, each distribution dP is a probability, i.e. 
is normalised to 1 when integrated on its domain of applicabil¬ 
ity. In terms of parametrisation: 


• as in the previous release, the spatial distribution 
d^y(r)/dy is selected from gENUM_PR0FILE (see table 1); 


• the mass distribution dPM(M)ldM is a simple power-law 
with two parameters (normalisation and slope ckm ~ 1.9), 
again similarly to the previous release; 


• the concentration distribution dPc(M,c)/dc is a new 
feature of the code (in clumps.h) and is chosen to 
be either a Dirac function or a log-normal distribution 
(gENUM_CVIR_DIST, see table 1): 


dc 


(M, c) : 


exp 


In c - ln(c(M)) 




^crc(M) 


C (Tc(M) 


( 10 ) 


^The concentration at a given characteristic overdensity A is defined to be 
CA = RAlr- 2 , where Ra is the radius of the clump for which the density equals 
this overdensity, and r _2 is the position where the logarithmic slope of the DM 
density profile of the clump reaches -2. 


U) 


r 

Jm. dM ^ X 


dM)dp 


Cmin(^) 


^ (M, c) £(M, c) dc dM, (13) 
dc 


and the average J value from the clumps becomes 


^ I 

Jo Jim 


(/subs) ~ Mot 


dPy 

dV 


(/, Q)d/dO 


JMmin 


dPM 

dM 


f 


x(M)d® 


Cmin 


dc 


(M,c)£(M,c) dcdM. 


(M) 


(14) 


The variance on J remains as calculated in Paper I, i.e. 


subs 




(15) 


but the difference is that {£^} and {£} must now be calculated 
accounting for the extra integration on the distribution of con¬ 
centration dPc(M,c)/dc. Using gENUM_CVIR_DIST=kL0GN0RM 
with a typical scatter of 0.18 around (Cyir) [20] leads to an in¬ 
crease of {£} by ~ 15% at all masses (see the documenta¬ 
tion in clumps. h for illustrative plots), compared to the use 
of gENUM_CVIR_DIST=kDIRAC (as in the previous release). 


3.1.3. Substructure boost: including extra-levels 

Based on the hierarchical structure formation scenario, the 
calculation of /subs and /cross-prod should include a multi-level 
description of the clumps, i.e. contributions from clumps within 
clumps within clumps, etc. 

To exemplify how this is implemented in CLUMPY, let us con¬ 
sider a DM halo, with density profile fully encompassed in 
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the integration angle. Not considering any substructure within 
this halo, its intrinsic luminosity is directly given by Eq. (12), 
denoted ‘0’ to highlight the absence of substructures in that DM 
halo: 


Xo(M,c) 


= r 

^ Vcl 


dV. 


(16) 


A hierarchy of n levels of substructures within this host halo 
{n= 1,2,3, etc.), may be computed from the n - I level as (we 
drop the implicit dependence on c for compactness): 


Xn(M) - Xsm(^) + Xcross-prod(^) (17) 


with 


+ 


NUM) 



(M) 


Xsm(M) 


X. 


'Cross-pro( 


,d(M) 


^Vcl 

2 r p:r(^)<Psubs(M))dy. 

-'Vel 


In the first version of the code, only n = \ was considered. 
This has now been extended to any level with a recursive im¬ 
plementation of Eq. (17). Part of this implementation requires 
to identify the radius at which the slope of the spatial distribu¬ 
tion of the clumps equals -2, i.e., it can only be enabled with 
distributions whose outer slope is steeper than -2. 

The impact of the substructures as well as the relative con¬ 
tributions of higher level subhaloes are illustrated in Eig. 1. 
The top panel corresponds to the boost factor X 4 /X 0 calculated 
from Eq. (17) as a function of the host halo mass Myir- The 
exact dependence depends on several key parameters (minimal 
mass of subhaloes, mass distribution slope au, and Cyk - ^vir 
parametrisation), but our results using, e.g., kSANCHEZ14_200 
for Cyir - /?yir are in agreement with the results of [22]. The 
bottom panel details the contributions of higher-order contribu¬ 
tions with respect to the previous level. The main contribution 
(not shown) is the first level of subhaloes (i.e. n = \), which 
is responsible for most of the boost seen in the top panel. The 
second level of substructures further contributes at most to 30% 
(solid red line), for the most massive haloes. As underlined in 
several previous studies using a different approach [23, 22], the 
third level of substructures contributes to less than 5% of the 
total. 


3.1.4. Recommendation regarding concentration and number 
of substructure levels 

Concentration parametrisations. The gDM_FLAG_CVIR_DIST 
parameter allows the user to choose the mass-concentration 
from eight pre-defined parametrisations taken from various 
studies [20, 22, 33, 34, 35, 36, 37, 38] and listed in Table 1. 
While any of these parametrisations may be used to charac¬ 
terised the subhaloes concentration in CLUMPY, it is important 
to keep in mind that, to date, these parametrisations have been 
established for main haloes only (i.e. not substructures). Eur- 
thermore, note that some of these parametrisations are sim¬ 
ple power-laws, shown to provide good fits to simulations in 
the mass range ~ 10^^ - lO^^M© [34, 35] or to X-ray galaxy 






Figure 1 : Top panel: boost factor X4 /Xo as a function of the host halo mass for 
two aM and two Cyir - ^vir relationships. Bottom panel: relative importance of 
higher-level substructure contributions. 


cluster data [36]; these should however not be extrapolated to 
lower masses where they would overestimate the concentra¬ 
tion. Parametrisations allowing for a flattening of the mass- 
concentration relation at lower mass such as [20, 22, 33, 38] 
are more realistic and the default (and recommended) CLUMPY 
setting is gDM_FLAG_CVIR_DIST=kSANCHEZ14_200 [22]. 

Number of substructure levels. The level parameter 
(gDM_SUBS_NUMBEROFLEVELS in clumpy_parains.txt) is 
set to ^ = 1 by default. Setting ^ = 0 is not implemented 
directly but is obtained instead by asking that the fraction of 
the mass of the host halo under the form of substructure is zero 
(gTYPE_SUBSJy[ASSFRACTION = Oinclumpy_params.txt). 

Note that enabling more than one level of substruc¬ 
tures (gDM_SUBS_NUMBEROFLEVELS in clumpy_params.txt), 
and/or requiring a log-normal distribution of concentrations 
(gENUM_CVIR_DIST=kL0GN0RM), significantly increases the 
computational time required to run CLUMPY. We find that us¬ 
ing the concentration scatter increases by ~ 30% the computa¬ 
tional time for a given sky map. Asking for a second level of 
substructures almost doubles the duration of the run. Requir- 


5 








ing simultaneously the concentration scatter and a second level 
of substructures increases by a factor ~ 10 the computational 
time. We underline that using n = 2 multi-levels is enough to 
reach a precision better than 5%. 

3.2. Host halo triaxiality 

Both numerical simulations and observations hint at triaxial 
rather than spherical DM haloes [49]. In this CLUMPY release, 
we supplement the default spherical halo configuration with a 
more general triaxial model (a b c), where a, b and c cor¬ 
respond to dimensionless major, intermediate and minor axes. 
Cosmological simulations have shown that these axes can vary 
as a function of the iso-density radius Ri^o [49], but this refine¬ 
ment is not considered here. The position (X, Y, Z) in the coordi¬ 
nate framework attached to a given DM halo centre corresponds 
to the iso-density radius 

y2 ^2 

In this configuration, the DM density is given by p(R), where 
p can be any of the spherical profiles already implemented in 
CLUMPY (e.g., Einasto, Zhao, ...). 

3.2.1. Different implementations for different DM halo types 

The treatment of triaxiality relies on seven parameters (three 
axes for the shape, three Euler rotation angles, for the orienta¬ 
tion, and one boolean to switch on or off triaxiality). 

• Galaxy: for the Milky Way, these parameters are denoted 
gGAL_TRIAXIAL_XXX and must be defined in the param¬ 
eter file clumpy_params.txt (see table 2). Whenever 
triaxiality in enabled, the options propagates to the total, 
smooth and average clump contributions. 

• Halo list: CLUMPY can also run on a user-defined 
list of haloes. In addition to the position and struc¬ 
tural parameters of each halo, already required in 
the first version of the code, the halo list may now 
include (as an option) the seven extra parameters 
for triaxiality (compare data/list_generic.txt and 
data/list_generic_triaxial. txt). Each halo in the 
list can therefore have its own triaxial properties. See the 
documentation for the required format and examples. 

• Substructures: triaxiality is not enabled for the description 
of subhaloes in its host halo, since it would require a sta¬ 
tistical distribution of axis ratios and orientations. Eirst, 
the orientation is not important if the integration angle en¬ 
compasses the subhalo volume (and this is the case for all 
subhaloes but a tiny fraction). Second, the exact distri¬ 
bution of axis ratios is not known, but the impact in the 
calculation is certainly sub-dominant compared with other 
uncertainties (concentration, profile, etc.). 


3.2.2. New CLUMPY functions 

Erom the line-of-sight integration point of view, no major 
changes were required as the first release was already designed 
to deal with non-spherical halo integrations (in the function 
los_integral 0), although this capability was not used at the 
time. Triaxiality only requires a few new functions: 

• get_riso_triaxial() to transforms (l,a,l3) integration 
position into (v,y,z) DM halo coordinates—and if re¬ 
quired (X, T, Z) Euler rotated target coordinates—from 
which 7?iso given in Eq. (18) is calculated; 

• mass_triaxialhalo() must be used instead of 
mass_singlehalo 0 to obtain the mass of a triaxial halo. 

These changes are transparent to the user, whose only concern 
must be to provide the seven triaxiality initialisation parame¬ 
ters. 


4. New outputs (content and format) and pixelisation 

4.1. y-ray and v spectra: spectra. h 

We add the gamma and neutrino spectra from dark matter 
annihilation and decay. We use the tabulated spectra of recent 
PYTHIA simulations [19, 50] including or not EW corrections 
[51], in which a Higgs mass of 125 GeV is assumed. 

4.1.1. Implementation of PPPC4DM1D in CLUMPY 

The values are calculated by a 2D linear interpolation on 
log(£’) and moM f^om tabulated spectra [51]. 

• Branching ratios: gPP_BRANCHINGRATIO_LIST (see ta¬ 
ble 2) is a list of comma-separated values for the branch¬ 
ing ratios^ of the 28 primary channels e^e^, e^e~, 

p"p~, TrTr, qq, cc, bb, tt, 

WjWj, W^W~, ZjZi, ZjZt, ZZ, gg, yy, hh, v^Ve, v^v^, 
yy ^ 4e, VV 4p, VV 4 t. 

• Einal state: gENUM_FINALSTATE (see table 1) must be 
chosen among kGAMMA (y-rays) and kNEUTRINO (v). 
The flavour of the latter is chosen (see table 2) among 
gSIMU_FLAG_NUFLAVOUR = kNUE, kNUMU, kNUTAU. 

• Decay or annihilation: it is set by the boolean 

gPP_DM_IS_ANNIHIL_OR_DECAY. For a decaying DM par¬ 
ticle, we use the same spectra dNy^yjdE as for the 
case of annihilation, but assuming = 

dN^’y'IdEim’oyilD. 


^For a minimum of safety, BRi < 1 is checked and indicated to the 
user but ZfJi BRi > 1 can be allowed due to the possible redundancy between 
channels with polarisations or chiralities. 
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Table 1: Enumerators and allowed keyword (and reference) in the CLUMPY code. Highlighted entries (shade of grey) correspond to new enumerators and/or keywords 
of this version (w.r.t. the first release). 


Enumerator 


Flags available 


gENUM_ANISOTROPYPROFILE kCONSTANT, kBAES [24], kOSIPKOV [25, 26] 

gENUM_LIGHTPROFILE kEXP2D [27], kEXPSD [27], kKING2D [28], kPLUMMER2D [29], kSERSIC2D [30], kZHA03D[31, 32] 

gENUM.CVIRMVIR kB01_VIR [20], kENS01_VIR [33], kNET007_200 [34], kDUFFY08F_{VIR, 200, MEAN) [35] 

kETT0RI10_200 [36], kPRADAll_200 [37], kGI0C0LI12_VIR [38], kSANCHEZ14_200 [22] 
gENUM_CVIR_DIST kLOGNORM [20], kDIRAC 

gENUM_PROFILE kZHAO [31, 32], kEINASTO [39], kEINASTO_N [40], kBURKERT [41], 

kEINASTOANTIBIASED_SUB [42, 43], kGA0_SUB [44, 45] 


gENUM_TYPEHALOES 
gENUM_FINALSTATE 
gENUM_NUFLAVOUR 
gENUM PP SPECTRUMMODEL 


kDSPH, kGALAXY, kCLUSTER 
kGAMMA, kNEUTRINO 
kNUE, kNUMU, kNUTAU 

(kBERGSTR0M98 [46], kTASITSI0MI02 [47], kBRINGMANN08 [48])^ , kCIRELLIll{EW, NOEW) [19] 


* Spectra for y-rays only. 


4.1.2. Neutrino mixing 

The simulations from [19] provide neutrino production spec¬ 
tra. For distant astrophysical sources, the journey in vacuum 
and transition between the different flavour states can be de¬ 
scribed by average oscillations [52]: 

3 

P{vi ^ w) = P(vi ^ n) = WnfWikf (19) 

k=l 

where U is the neutrino mixing matrix and k = 1,2,3 

for the 3 mass eigenstates. The oscillation matrix is filled 
with nu_oscillationmatrix() from the mixing angles 
{0i2, 6 > 23 , Ou) (gPP_NUMIXING_THETA{12,13,32}_DEG (see ta¬ 
ble 2), whose default values are taken from [53], i.e. 
{34°, 49°, 9°). The CP phase is taken to be zero (see [53] for 
a recent discussion). Oscillation effects are applied to the v^, 
and Vj spectra. The code gives the resulting v^, or Vj fiuxes^^. 

4.1.3. Particle physics term and flux in CLUMPY 

The above final states, the particle physics term Eq. (2), 
or the y-ray (or neutrino) flux (Eq. 1) for a given astrophys¬ 
ical factor can be displayed with ./bin/clumpy -z. The 
corresponding functions in CLUMPY are dNdE(), dPPdEO, 
dPPdE_integrated(), and flux(). Fluxes (and integrated 
fluxes) are also displayed whenever the astrophysical fac¬ 
tor for the Galaxy (./bin/clumpy -g) or for a DM halo 
(. /bin/clumpy -h) is calculated. 

4.2. Map handling and outputs: healpix_f its .h 

This new release provides additional tools in the context of 
2D maps, as well as more advanced output options in addition 
to the ASCII files outputs and plots which were available with 
the first version. 


^^Relevant quantities for DM detection with neutrino telescopes are + 
y^ = 2 X y^. This factor 2 for anti-neutrino contribution is not accounted for in 
CLUMPY results. 


4.2.1. YlEkCPlyi map tools 

CLUMPY is now interfaced with the HEALPix (Hierarchical 
Equal Area isoLatitude Pixelation) package, which provides 
a large set of routines for efficient manipulation and analysis 
of numerical data on the discretized sphere [54]. Within the 
HEALPix framework, CLUMPY now robustly handles large field 
of views (FOV) at arbitrary positions of the sky, and the fast 
calculation of full skymaps. 

Various shapes for the part-sky grid are now available 
(i.e., circular, great-circle or iso-latitude/longitude rectangular 
shapes of the FOV). Reversely, for full-sky calculations, sim¬ 
ple masking shapes can be applied to reduce computation time, 
memory usage and output file size of the simulation. Examples 
for available FOV choices and masks are given in the Doxygen 
documentation. 

Using the HEALPix library routines, CLUMPY now provides 
the options of smoothing the output /-factor skymaps with a 
Gaussian beam and calculating the angular power spectrum 
(APS) of the maps. Up to two smoothing kernels may be spec¬ 
ified (e.g., one for a y-ray instrument and one of a neutrino 
observatory), both being applied on the same output map origi¬ 
nally computed by CLUMPY^The APS calculation^^ is limited 
to the Galactic mode (. /bin/clumpy -g) and is performed in¬ 
dependently for the total /-factor profile of the halo /tot, the 
smooth halo contribution /sm, the substructure component /subs 
(average (/subs) and drawn haloes /drawn contributions) and the 
cross-product /cross-prod (annihilations only)^^. 

4.2.2. Output format 

The output files of the 2D skymaps and APS are now written 
in the fits (Flexible Image Transport System) file format [55]. 


order to avoid smoothing the edges of part-sky or masked maps, CLUMPY 
appropriately extends the original grid when smoothing is requested by the user. 

^^Note that power spectra for part-sky and masked maps are affected by 
power suppression and spectral leakage effects according to the shape of the 
windowing function of the FOV. Further accounting and correction for these 
effects must be performed outside of CLUMPY. 

^^Note that the power spectrum Q of the /-factor components is an additive 
quantity, Cf^ + ... . 
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Table 2: CLUMPY parameters for the user-defined input parameter file (clumpy_params. txt). For the sake of completeness, all parameters are reproduced sorted 
by block (cosmology, dark matter, particle physics, etc.), with new parameters highlighted in grey, and deprecated ones in strikethrough. 


Name 


Definition 


Cosmological parameters (updated from Planck results) 


gC0SM0_HUBBLE 

gC0SM0_RH00_C 

gC0SM0_0MEGA0_M 

gC0SM0_0MEGA0_LAMBDA 


Hubble expansion rate h - km s"^ Mpc"^) [-] 

Critical density of the universe [M© kpc"^] 

Present-day pressure-less matter density 
Present-day dark energy density 


Dark matter parameters 

gDM_FLAG_CVIR_DIST 

gDM_LOGCVIR_STDDEV 

gDM_SUBS_NUMBEROFLEVELS 

gDM_MMIN_SUBS 

gDM_MMAXFRAC_SUBS 

gDM_RHOSAT 


Distribution around c{M) from which concentrations are drawn: {kLOGNORM, kDIRAC} 
Width of log-normal c{M) distribution (if gDM_FLAG_CVIR_DIST=kLOGNORM) 

Number of levels for subhaloes 
Minimal mass of DM haloes [M©] 

Defines the maximal mass of clump in host halo: M^ax = gDMjy[MAXFRAC_SUBS X Mhost 
Saturation density for DM [M© kpc"^] 


Generic (sub-)halo structural parameters (TYPE = DSPH, GALAXY or CLUSTER) 

gTYPE_CLUMPS_{FLAG_PROFILE, . . .} Description of subhaloes for host TYPE: c(M), inner profile, shape parameters 

gTYPE_DPDM_SLOPE Slope of the clump mass function 

gTYPE_DPDV_{FLAG_PROFILE, RSCALE, . . .} Spatial distribution of substructures in object TYPE 
gTYPE_SUBS jyiASSFRACTION Mass fraction of the host halo in clumps 


Milky-Way DM (sub-)halo structural parameters 


gGAL_CLUMPS_{FLAG_PROFILE, 
gGAL_DPDM_SLOPE 

gGAL_DPDV_{FLAG_PROFILE, RSCALE, 
gGAL_SUBS_{Ml, M2, N_INM1M2} 
gGAL_{RH0S0L, RSOL, RVIR} 
gGAL_TOT_{FLAG_PROFILE, RSCALE, . 
gGAL_TRIAXIAL_AXES [0-3] 
gGAL_TRIAXIAL_ROTANGLES [0-3] 
gGAL_TRIAXIAL_IS 


Description of Milky-way DM subhaloes 
Slope of clump mass function 
Spatial distribution of substructures in object TYPE 
Number of Milky-Way subhaloes in [Mi, M 2 ] 

Local DM density [GeV cm"^], distance GC-Sun [kpc], virial radius [kpc] 
Description of the total DM profile 

Dimensionless major (a), intermediate (Z?), and minor (c) axes (see Eq. (18)) 
Euler rotation angles for triaxial Milky-Way halo [deg] 

Switch-on or off triaxiality calculation (i.e., use or not the 2 parameters above) 


Particle physics iugredieuts (for y-ray and v flux calculation) 

gPP_BR [gN_PP_BR] List of comma-separated values of branching ratios for the 28 channels 

gPP_DM_ANNIHIL_DELTA Eor annihilating DM, factor 2 in calculation if Majorana, 4 if Dirac 

gPP_DM_ANNIHIL_SIGMAV_CM3PERS Eor annihilating DM, velocity averaged cross-section {crv)o [cm^ s"^] 

gPP_DM_DECAY_LIFETIME_S Eor decaying DM, lifetime Tdm of DM candidate [s] 

gPP_DM_IS_ANNIHIL_OR_DECAY Switch for annihilating or decaying DM {replace deprecated gSIMU-IS-ANNIHIL_OR_DECAY) 

gPP_DMjy[ASS_GEV Mass moM of the DM candidate [GeV] 

gPP_FLAG_SPECTRUMMODEL Model to calculate final state {replace deprecated gDM_GAMMARAY_FL AG-SPECTRUM) 

gPP_NUMIXING_THETA{12, 13, 23}_DEG Neutrino mixing angles [deg] 

Simulation parameters/outputs (for a given CLUMPY run) 

gLIST_HAL0ES DM haloes considered in /-factor calculations [default=dat a/list .generic. txt] 

gLIST_HALOES_JEANS Objects considered in Jeans’s analysis [default=dat a/list .generic, j eans. txt] 

gSIMU_ALPHAINT J)EG Integration angle aint [deg] (if gSIMU.HEALPIX.NSIDE not -1, use HEALPix resolution) 

gSIMU_EPS Precision used for any operation requiring one (numerical integration, ...) 

gSIMU_SEED Seed of random number generator to draw clumps (if 0, from computer clock) 

gSIMU_FLAG_NUFLAVOUR Choice of neutrino flavour (kNUE, kNUMU, kNUTAU) 

gSIMU_FLUX_AT_E_GEV Energy (GeV) at which to calculate fluxes 

gSIMU_FLUX_E.MIN Lower energy bound (GeV) for the integrated flux calculation 

gSIMU_FLUX_E.MAX Upper energy bound (GeV) for the integrated flux calculation 

gSIMU_GAUSSBEAM_GAMMA.FWHM.DEG Gaussian beam [deg] for y-ray detector for skymaps smoothing (no smoothing if set to -1) 

gSIMU_GAUSSBEAM_NEUTRINO_FWHM_DEG Gaussian beam [deg] for y detector for skymaps smoothing (no smoothing if set to -1) 

gSIMUJlEALPIX_NSIDE Aside of HEALPix maps (if -1, set to be as close as possible to aint) 

gSIMUJlEALPIX_RING_WEIGHTS_DIR Ring weights directory for improved quadrature (optional) 

gSIMU_IS_ASTR0_0R_PP_UNITS Outputs (plots and files) in astro (M© and kpc) or particle physics (GeV and cm) units. 

gSIMU_IS_WRITE_FLUXMAPS Eor 2D skymaps, whether to save or not y-ray and v fluxes (the / factor is always saved) 

gSIMU_IS_WRITE_FLUXMAPS_INTEG_ORJ)IFF If gSIMU_IS_WRITE_FLUXMAPS is true, whether to save integrated or differential fluxes 
gSIMU.IS.WRITE.GALPOWERSPECTRUM Whether to calculate (and save) or not the DM power-spectrum for the Milky-Way 

gSIMU-IS.WRITEJlOOTFILES Whether to save or not .root files even if option -p is used (not enabled for skymaps and ’stat’) 

gSIMU OUTPUT DIR Output directory to select other than local run (directory is output / if set to -1) 
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This format allows to clearly and efficiently store several maps 
as binary tables in one single file, and to append extensive meta 
information in the fits headers of the output file. Both full- 
sky and part-sky HEALPix maps can be handled in the fits 
format. CLUMPY fits files are composed of several extensions 
(/-factor skymaps, smoothed /-factor skymaps, fluxes), each 
containing different contributions (/sm, /subs, ^cross-prod.- - -) in 
separate columns of the extension. Each extension possesses 
its own fits header, containing the relevant input parameters 
of the performed simulation. Writing fits file in CLUMPY is 
done by using the cf itsio library either directly, or indirectly 
via existing HEALPix routines. 

Many programs are available for reading, displaying and 
manipulating HEALPix maps in the fits format^"^ (e.g. fv^^, 
DS9^^, Aladin^^, Healpy from the HEALPix package, ...). To 
facilitate the display and post-processing of 2D-simulations, the 
. /bin/clumpy -o mode has been implemented to convert and 
extract the output maps in different file formats (more details 
are given in the README CLUMPY file). In particular, ASCII 
files similar to those of the first CLUMPY release can be writ¬ 
ten. For 2D ROOT plots and output, the data is transformed 
from the HEALPix pixelisation scheme onto a rectangular grid 
in longitudinal and latitude coordinates. This tangential plane 
approximation is only appropriate for small FOV. 

5. New science: Jeans analysis 

The original goal of the CLUMPY code was to provide a robust 
and versatile tool to compute the astrophysical factors for DM 
annihilation and decay signals. The updates and new features 
introduced to improve CLUMPY in that respect were presented in 
§3. Here, we move to an entirely new aspect of this release, 
namely the possibility to run a Jeans analysis coupled to an 
MCMC engine, in order to derive data-driven DM density pro¬ 
files from a set of stellar/galactic kinematic data. This method 
is independent from the /-factor calculation itself. The Jeans 
module presented here has already been used in [56], [57] and 
[58], in order to provide robust /-factors (and associated error 
bars) for both the ‘classical’ and ‘ultrafaint’ dwarf spheroidal 
galaxies of the Milky Way, which are among the most promis¬ 
ing targets for DM indirect detection. The recently discovered 
dSph galaxy candidates [59, 60] are typical objects on which to 
apply the Jeans analysis once kinematics data become available, 
in order to obtain their /-factor. 

5.7. Jeans equations 

The Jeans equation is obtained after integrating the colli- 
sionless Boltzmann equation in spherical symmetry, assuming 
steady-state and negligible rotational support [18]. It relates the 
dynamics of a collisionless tracer population (e.g. stars in a 
dwarf spheroidal galaxy or galaxies in a galaxy cluster) to the 


^^http://fits.gsf c.nasa.gov/fit s_viewer.html 
^^http://heasarc.gsfc.nasa.gov/ftools/fv/ 
^^http://ds9.si.edu 
^^http://aladin.u-strasbg.fr 


underlying gravitational potential. Namely, the Jeans equation 
reads 


^ d I GM(r) 

V dr ^ ' r 


( 20 ) 


with 


M(r) = An 


f 


Ptot(s)s^ds , 


( 21 ) 


and where: 


• the definition of the enclosed mass M(r) assumes that the 
tracer population (i.e. the stars) has a negligible mass com¬ 
pared to the underlying DM halo. The density profile of 
the latter is Ptot(^)’, 


• y(r) is the 3D number density (or light profile) of the tracer 
population. Its 2D projection is called surface brightness 
(see below); 

• is the radial velocity dispersion of the tracers; 

• y^ani(^) = 1 - v^/v^ is the velocity anisotropy of the trac¬ 
ers, and depends on the ratio of the tangential to the radial 
velocity dispersions. 


The formal solution to the 3D Jeans equation is 

v(r)v2(r) = X r f(s)v(s)^^^^ds , (22) 

f(r) Jr 


with 


/(r) = fr^ exp 





(23) 


In practice, observations provide only the 2D projections on the 
sky of the tracer velocity dispersion and number density. For a 
given projected radius R, the projected 2D solution of the Jeans 
equation is 


mcr^piR) 


’f 


1 -/3am(i')^\y(r)v^r(0dr , 


(24) 


with I(R) and crp(R) the surface brightness and projected veloc¬ 
ity dispersion, respectively. A review on recent developments 
and applications of the Jeans analysis can be found, e.g., in [61]. 


5.2. Implementation in the new library j eans_analysis. h 

Computing Eq. (24) requires three levels of imbricated in¬ 
tegrations which can be time consuming. For some specific 
parametrisations of ySani (05 analytical shortcuts may be found 
to reduce the problem to a single integration. These shortcuts, 
based on the kernel implementation of Eq. (24) of [62], are im¬ 
plemented in CLUMPY, along with the full numerical integration 
whenever necessary. This approach is fully described in the 
documentation attached to the j eans_analysis. h library and 
is not detailed again in this paper. 

The functions at the core of the library are related to the de¬ 
scription of the ingredients of the Jeans analysis: 
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• beta_anisotropy(gENUM_ANISOTROPYPROFILE, . . .): 
returns the velocity anisotropy profile selected by the 
user in the file containing the halo structural parameter 
description. The most generic form for the anisotropy is 
the Baes & Van Hese parametrisation [24], 


y^ani(r) = 


ygp +ygoo(r/rj^ 
1 + {rlva)^ 


(25) 


with four free parameters (J3o,l3oo, ra and rf). The user may 
also choose the constant parametrisation or the Osipkov- 
Meritt parametrisation [25, 26], which both are special 
limiting cases of the Baes & Van Hese profile (see Ta¬ 
ble 1). 


• light_profile(gENUM_LIGHTPROFILE, . . .): returns 
the 2D or 3D light profile (i.e. 1{R) or y(r)) selected by the 
user in the file containing the light profile structural param¬ 
eters. Four 2D profiles and two 3D profiles are available 
(see Table 1). Depending whether the user selects a 2D or 
a 3D light profile, the corresponding (de-)projections are 
performed by CLUMPY in order to solve the Jeans equation. 

• jeans_*: a series of functions which provide the neces¬ 
sary steps to effectively solve the Jeans equation, through 
the various integrations detailed in §5.1. 


5.3. Jeans analysis within CLUMPY.* inputs and outputs 

The standard procedure of a Jeans analysis is the following: 

1. fit the surface brightness profile I{R) with a given paramet¬ 
ric function^ 

2. choose a parametrisation for the DM density profile ptot(^) 
and the tracer anisotropy I3^^i{r)\ 

3. compute the projected velocity dispersion o-p from I{R), 
Ptot(r) andy5ani(r) (Eq. 24); 

4. find the DM density and anisotropy parameters that best 
reproduce the kinematic observables. 

The DM density profile can then be used to compute the J- 
factor, or any other DM-related quantity. Figure 2 summarises 
the steps of the Jeans analysis, and the related executables. We 
describe below the different ingredients of the analysis. 


5.3.1. Input files: kinematic data and free parameters 
Kinematic data. The kinematic data of the tracer population 
(i.e., stars in a dwarf spheroidal galaxy) are usually in the form 
of velocity dispersion profiles (Tobs (resp. squared velocity dis¬ 
persion, or line-of-sight velocities v/. These data must 

be filled in specific data files. Three examples are given in the 
data/ directory, for the three different types of kinematic data: 
data_sigmap.txt, data_sigmap2.txt, and data_vel.txt. 
A keyword is associated to each data type and written in the 
header of the file: Sigmap, Sigmap2, and Vel. Note that sur¬ 
face brightness profiles can also be read by CLUMPY (see the 
example file data_light. txt). 


Free parameters. The surface brightness, velocity anisotropy 
and DM density profiles used in the Jeans analysis are set in 
the parameter file params_jeans.txt, located in the /data 
directory. Up to nine parameters can be let free in the analysis 
to describe velocity anisotropy and DM density profiles. Their 
ranges are set in params_j eans. txt. The parameters describ¬ 
ing the surface brightness profile, as well as some characteris¬ 
tics of the object (e.g. name, distance,...), must also be filled in 
this parameter file. 


5.3.2. Likelihood functions 

In j eans .analysis . h, we provide two likelihood functions 
log_likelihood_jeans_{binned,unbinned}0 for the fit of 
the velocity dispersion o-p to the kinematic observables. If the 
kinematic data are binned velocity dispersion profiles, the like¬ 
lihood function is^^: 


£ 


bin 


1^1 Acr,(/?,) 21 ) / 


(26) 


with IscrfRi) the error on the velocity dispersion at the radius Rt. 
If the data are in the form of line-of-sight individual velocities 
V/, then the likelihood function used is [63]: 


£ 


unbin 


^ ^ tracer' 

n 




(2ny 


- 1/2 


+ A2 


exp 


1/ (vj-v)^ V 
2Vo-p(^,) + A2 / 


(27) 


with Pi the membership probability of the i-th tracer^®. The 
underlying assumption is that the line-of-sight velocity distri¬ 
bution is Gaussian, centred on the mean velocity of the tracers 
V, with a dispersion of velocities at radius Ri of the i-th tracer 
coming from both the intrinsic dispersion crp(Ri) from Eq. (24) 
and the velocity measurement uncertainty Ap.. 


5.3.3. MCMC and bootstrap analyses 

With CLUMPY, the Jeans analysis can be done in two ways. 
The first uses the MCMC toolkit GreAT, and the second a sim¬ 
ple minimisation with a bootstrap analysis. 

MCMC analysis. The MCMC technique, based on Bayesian 
parameter inference, allows for an efficient sampling of the 
posterior probability density function (PDF) of a vector of pa¬ 
rameters with Markov chains. Credibility intervals (CIs) can 
then be reconstructed from the PDFs, for any quantity of inter¬ 
est. CLUMPY is interfaced with the Grenoble Analysis Toolkit 
(GreAT) [65], a public C++ statistical framework which han¬ 
dles MCMC analysis with the Metropolis-Hastings algorithm 
[66]. To use GreAT with CLUMPY, the user has to create the en¬ 
vironment variable ‘GreAT’, pointing to the installation folder 
of GreAT. The jeansMCMC executable provides an example of 
MCMC analysis, using the input files and likelihood functions 
described in the previous sections. 


While CLUMPY does not perform the actual fit to the surface brightness pro¬ 
file, several parametrisations of I(R) are implemented in the code as to perform 
the Jeans analysis (see gENUM_LIGHTPROFILE enumerator in Table 1). 


^^If the kinematic data are binned squared velocity dispersion profiles, then 
cr is replaced by cr^ in Eq. (26) 

Pi estimates the probability that the i-th. tracer belongs to the object. This 
quantity is often available for stars in dwarf spheroidal galaxies. 
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Bootstrap analysis. The Jeans analysis can also be run with a 
simpleminimisation, using the Minimizer class of the ROOT 
libraries. The;^^ function is 

X^ = -2\ogU), (28) 

with X being one of the likelihood functions defined in the pre¬ 
vious section (Eqs. (26) or (27)). The jeansChi2 executable 
provides an example of;^^ analysis. 

With the same executable, the user can also run a boot¬ 
strap analysis for an estimation of the statistical error bars [67], 
which proceeds as follows: i) for a given kinematic data sam¬ 
ple of n line-of-sight velocities, N re-samples are generated by 
drawing with replacement n velocities among the original sam¬ 
ple; ii) for each re-sample, the best-fitting anisotropy and DM 
parameters are found using the minimisation; iii) the boot¬ 
strap estimator of any quantity of interest is the mean over the 
N configurations, and the dispersion is used as statistical uncer¬ 
tainty. 


data_sigmap. txt, data_sigmap2. txt, 
data vel.txt 


Jeans data 


pa ramsJeans.txt 


Jeans 

parameters 


Input files 


j eans_analysis.h 



I Jeans module 


stat_example.dat 


'Statistical' 

file 


./clumpy -s 


5. 3.4. 'Statistical ’ output file 

Each Jeans executable creates an output ‘statistical’ file, 
readable by CLUMPY with the option -s. The latter allows to 
draw the estimators of several DM-related quantities, such as 
the astrophysical factors or the density profile, from the results 
of the statistical analysis (i.e., median and credible intervals 
from an MCMC analysis, mean and and dispersion for a boot¬ 
strap analysis). With respect to the previous CLUMPY version, 
new options were added to draw several Jeans-related quanti¬ 
ties, i.e. o-p{R), I{R), etc. 




-N 

J,D 

V_ ) 


PdM'>^ 


f -\ 

\ _/ 


Output files and ROOT 
plots 


Figure 2: Diagram of the Jeans analysis with CLUMPY. From a kinematic data 
file and a parameter file describing the free parameters, a statistical Jeans analy¬ 
sis can be run with . / j eansMCMC or . / j eansChi2. A ‘statistical’ output hie is 
created, from which estimators of different quantities of interest (i.e., /-factors) 
can be obtained with . / clumpy -s. 


6. Installation and run examples 

In this section, we provide a quick description of CLUMPY’s 
architecture and installation, the various options for running the 
code and, finally, a few run and plot examples obtained with 
CLUMPY. Many more examples and illustrations are provided 
along with the Doxy gen documentation. 

6.1. Code installation!compilationiexecutables 

CLUMPY is written in C/C-f-f and relies of the ROOT, cf it sio, 
HEALPix (C-f-l- and E90), and GSL libraries. Instructions for in¬ 
stallations and links to downloads are given in the README file. 
The code’s structure is standard, with separate directories for 
the various pieces of code: declarations are in include/*.h, 
sources in src/*.cc, compiled libraries, objects and executa¬ 
bles are respectively in the lib/, obj/, and bin/ directories. 
With respect to the first version, two new libraries were cre¬ 
ated (jeans_analysis .h and healpix_f its .h), one greatly 
expanded (spectra. h), and two new executables added for the 
Jeans analysis (j eansChi2. cc and j eansMCMC. cc). 

The various executables (and options) available for a Jeans or 
/-factor related CLUMPY run should be self-explaining (just type 
the command lines below). See the next section for examples. 


6.1.1. Jeans analysis 

As detailed above (5.2), the Jeans analysis allows to extract 
the DM density profile of a DM halo from the kinematic data 
of its tracer population. 

• ./bin/jeansChi2: Jeans analysis with simple min¬ 
imisation and bootstrap analysis. 

• ./bin/j eansMCMC: Jeans/MCMC analysis (only if 
GreAT is installed). 

6.1.2. J-factor (and flux) analysis 

The list of available options have not changed much, but they 
include all the refinements described in this paper (for the cal¬ 
culation, outputs, and displays). 

• ./bin/clumpy -g[option]: /-factor in Galaxy for 
smooth and clumps (mean or drawn). 

• ./bin/clumpy -h[option]: /-factor in any halo for 
smooth and clumps (mean or drawn). 

• ./bin/clumpy -o[option]: CLUMPY skymap file ma¬ 
nipulation. 

• . /bin/clumpy -s [option]: PDE and confidence levels 
on p(r), /, CTpiR)... 
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• ./bin/clumpy -z: spectra, particle physics term, and 
flux (annihilation, decay). 

• ./bin/clumpy -f: append extension with flux maps to 
existing FITS tile. 

6.2. Sky patches, smoothing, y-rays and v 

For illustration of the new features available for the 2D- 
skymap mode (-g5, -g8 and -h4, -h5) of CLUMPY, we pro¬ 
vide example plots for Milky-Way like DM haloes in Figs. 3, 4 
and 6. In these examples, the total density profile of the halo 
is parametrised by an Einasto profile with r _2 = 15.14kpc, 
aE = 0.17 [70] and a local dark matter density of p© = 
0.4 GeV cm"^ at R© = 8.0kpc. For a typical demanded preci¬ 
sion of RE < 2% (relative error), clumps down to masses 
of ~ 1 M© are drawn depending on the resolution, skymap size, 
... (see Paper I [10] and especially Eq. (28) in there for further 
details). The spatial distribution of subhaloes d^y/dV is de¬ 
scribed by the kGA0_SUB profile, 

N{<x) {\+ac)xP r 

- = —-, V = - , (29) 

A^200 \ + acx^ r2oo 

with ac = 11, a = 2, JS = 3 and r 2 oo = Rvir = 260 kpc. 
Note that for such values, the outer slope of the subhaloes spa¬ 
tial distribution is < 2, which does not allow the inclusion of 
more than one level of substructures as described in Sect. 3.1.3. 
The mass distribution &Pm/6M follows a power-law with in¬ 
dex aM = 1.9, normalised by an abundance of 150 subhaloes in 
the mass range [10^ M©, 10^^ M©]. The Cyk - Rvir relationship 
is modelled by kSANCHEZ14_200, with a Gaussian log-norm 
scatter of = 0.14 around (Cyir). A boost of the smooth J- 
factor is calculated taking into account substructures down to a 
minimal mass of 10“^ M©. Only one substructure level is ac¬ 
counted for in these examples. For a more detailed description 
and listing of the input parameters entering these simulations, 
we refer to the documentation and the standard parameter file 
clumpy_parains. txt delivered with the code. 

Figure 3 presents full skymaps for an observer located at 
R© = 8.0kpc, looking towards the centre of the galactic halo 
(-g7 mode). The maps are given in terms of differential y-ray 
intensities at 4 GeV for annihilation into XX channel and 
= 200 GeV. We provide an example for a spherical and a 
triaxial host halo shape. Flux and intensity maps both for y-ray 
and neutrinos can be generated from the the /-factor maps via 
the -f option, or directly together with the /-factor calculation 
when specified in the input parameters. 

Figure 4 shows the spherical halo within a small patch of the 
sky computed with the highest possible resolution available in 
the HEALPix-scheme, Vside = 2^^, together with the same patch 
smoothed with an instrumental beam by the HEALPix smooth¬ 
ing algorithm implemented in CLUMPY. Here, the /-values are 
given per steradian; /-factors for arbitrary opening angles equal 
or larger than the grid resolution can be derived from this by 
integrating over the opening angle AQ. For such high resolu¬ 
tions, the smoothing via spherical harmonics - as done by the 
implementation in CLUMPY - is very computationally expensive 



1 X 10-^0 1 X 10“^ 1 X 10“^ 1 X 10“^ 

[cm“^ sr“^ GeV“^] 


Figure 3: Top: Differential intensity skymap of the full galactic halo for y-rays 
from annihilation at 4 GeV for = 200 GeV and xX channel. The 

skymap is drawn in the -g7 mode for a numeric resolution of Vside = 2^ (cor¬ 
responding to a pixel diameter of ~ 0.12°) with the parameters for a spherically 
symmetric halo specified in the text. For a relative error RE< 2% (see 
Eq. 28 in [10]) and a numeric precision of 1%, a total number of 164,186 sub¬ 
structures is drawn within 7 h, consuming 1.0 GB RAM. Bottom: Galactic halo 
computed with the same parameters as on the top panel, but with a triaxially 
distorted shape of the total halo. The semi-axes ratios are, motivated by [69], 
set to bla = 0.83, c/a = 0.67, and a = 1.47. The spheroid is rotated versus the 
line of sight by the Euler angles a = 30°, p = 40° and y = 20° (arbitrary val¬ 
ues chosen for illustration only). This skymap needs 3/? of computation time 
(820 RAM) for 36,994 drawn substructures. However, as the spherical 
symmetry of the halo is now broken, a substantial amount of computation time 
(~ 1 hr) is required for the computation of the smooth /-factor components. 

(~ hours with > 20 GB RAM). For small FOV, smoothing via 
2-dimensional Fourier transform in the tangential plane approx¬ 
imation is recommended instead; this feature is not yet imple¬ 
mented in CLUMPY. 

Figure 5 shows two examples of plots generated by CLUMPY 
when populations of DM haloes are available (e.g. from a user 
list or for clumps drawn in the Galaxy). For more illustrations, 
see the population study performed with CLUMPY on clusters of 
galaxies in [15, 16, 17] and in CLUMPY runs. Note that the list 
of subhaloes is stored in ASCII flies (extension . drawn). 

To better compare the results of CLUMPY with results directly 
derived from the Aquarius [68, 70, 71], Via Lactea II [72, 43] 
and similar simulations of galactic haloes [73], Figure 6 dis¬ 
plays the intensity angular power spectrum (APS) of the signal 
over the full sphere. The angular power spectrum Q of an in- 
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Figure 6: Intensity APS for galactic substructures on the full sky for two differ¬ 
ent concentration models. For the particle physics term, y-rays from annihila¬ 
tion at 4 GeV for = 200 GeV and^^ ^ bb channel are assumed. The lines 
denote the mean value for 100 realisations of the same halo, the bands give 
the corresponding Icr containment range. The spectrum of the smooth halo 
7sm (grey dotted-dashed line), and the average contribution of the substructures 
((•^subs) + -^cross-prod5 solid green line), are plotted for comparison. 


Figure 4: Top panel: Example of a 5° x 5° skymap towards the galactic anti¬ 
centre of the spherically symmetric halo, obtained from the skymap mode 
-g7 for a numeric resolution of Aside = 2^^ (corresponds to a pixel diame¬ 
ter of ~ 0.007°). The colour scale gives the /-values per steradian in case 
of annihilation. This skymap contains 64,763 drawn substructures, computed 
with both a numeric precision and a relative error REj , < 1 % within 2,5 

hours (600 MB RAM), using the physical parameters specified above. Bottom 
panel: Same skymap as on the top, but smoothed with a Gaussian beam of 
FWHM = 0.1°. 



Figure 5: Example of two population study plots obtained for galactic sub¬ 
haloes towards the anti-centre in a region of size 45° x 45° (./bin/clumpy 
-g7 clumpy_params.txt 180. 0. 45. 45. 1). Left panel: Num¬ 

ber of substructures drawn (above the DM continuum) as a function of their 
apparent size as = (r^ is the scale radius, and d the distance from the ob¬ 

server). Right panel: cumulative of the signal of drawn subhaloes above a given 
/ (blue), and corresponding number n of subhaloes contributing to this signal. 
The red curve corresponds to the cumulative from n background regions. 


tensity map /(t^, (f) is defined as 

m 

with aim the coefficients of the intensity map decomposed into 
spherical harmonics 7/^, 

^max fn=+£ 

(f) = zz (31) 

£=0 m=-t 

The filled bands correspond to the power spectrum of 
the drawn substructures, when considering a spherical host 
halo (using the same parameters as in Fig. 3) but two 
different parametrisations of the Cyk - Rviv relationship 
(kSANCHEZ14_200 and kB01_VIR). Each spectrum was com¬ 
puted for 100 different statistic realisations of the same galactic 
halo. The bands in Figure 6 indicate the Icr containment ranges 
around the mean (assuming a log-normal variation of the Q) 
of the cosmic variance of the simulated galactic halo. The total 
average contribution of the substructures ((/subs) + /cross-prod. 
Eqs. (14) and (7)) is given by the green solid line. Apart 
from the largest angular scales, the drawn substructures com¬ 
pletely dominate over the APS of the averaged substructure 
contribution. The smooth signal (dotted-dashed line) dom¬ 
inates everywhere but for the smallest angular scales. The 
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Figure 7: Correlation plot obtained with the option -s2 of CLUMPY. An MCMC 
analysis was run previously with the default options of the jeansMCMC exe¬ 
cutable. 


kSANCHEZ14_200 model agrees well with the results presented 
in [70], where the resolved subhaloes from the Aquarius A-1 
simulation [68] are taken into account. These authors estimate 
the error on the APS, due to not resolving haloes below 10^ Mq, 
to be < 10%. In our approach, the full mass range for the sub¬ 
haloes (down to 10“^ Mo) is taken into account, within the /crit 
criterion, to ensure the required user-defined accuracy on the 
total /-factor. Testing various values of the latter, we find that a 
10 % accuracy requirement on the total /-factor translates into 
a ~ 1% error on the APS. 

6.3. Jeans: MCMC, PDF, CIsforp, J, cr^ 

Many useful quantities can be computed with CLUMPY after a 
Jeans analysis, using the statistical options (-s) of the CLUMPY 
executable: 

• Figure 7 shows a correlation plot made with the option 
-s2. These results were obtained with the default options 
of the jeansMCMC executable, i.e. using an MCMC anal¬ 
ysis with the GreAT package. The z-axis of the 2D plots 
shows the value of log £. 

• Figure 8 displays two examples of quantities obtained with 
clumpy -s after a Jeans analysis. The figure shows the 
median value and CIs on crp(R) and Jia-m), obtained with 
the default options of clumpy -s8 and clumpy -s6 re¬ 
spectively. The statistical output file used was obtained 
with the default options of the jeansMCMC executable and 
uses data of a mock ultra-faint dSph galaxy as provided 
with CLUMPY. For outputs obtained on real data, we refer 
the reader to Figs. 1 and 2 of [57]^^. 


^^The Jeans analysis implemented in CLUMPY is based on a data-driven 




Figure 8: Two examples of results obtained with the Jeans/MCMC analysis. 
The figure shows the median value and CIs on crp{R) and /(uint), obtained with 
the default options of clumpy -s8 and clumpy -s6 respectively. 


7. Conclusions 

The first version of the CLUMPY code [10] was originally de¬ 
veloped to provide a fast, robust and versatile tool to compute 
/ and D-factors for any halo, in a variety of configurations, and 
with up to one level of substructures. In this second release, sev¬ 
eral physically-motivated extensions have been added, such as 
drawing halo concentrations from a distribution around an av¬ 
erage mass-concentration relation, including triaxiality for DM 
haloes, allowing for several levels of substructures in the boost 
calculation, and providing new analysis tools (e.g., /-factor 
“population” studies). Furthermore, the code now includes the 
PPPC4DMID calculation for y-ray and neutrinos yields, allow¬ 
ing the computation of actual y-ray and neutrino differential or 
integrated fluxes. Finally, the 2D mode (i.e. skymaps) is now 
handled in the HEALPix pixelisation scheme for an improved 
behaviour over large fractions of the sky. 

Along with these extensions, a new module was also devel¬ 
oped to perform Jeans analyses on a set of kinematic data (e.g. 
from dwarf spheroidal galaxies). This allows the user to con¬ 
trol the entire analysis chain, from the reconstruction of the DM 


approach where as little prior from simulations as possible is used. In [57], we 
have compared our results to that of the Fermi collaboration and did not find 
any systematic shift between the /-factors of the two analyses. However, our 
data-driven approach results in larger CIs for the ultra-faint dSph galaxies given 
the very few kinematic data available for these objects (see figure 6 of [57]). 
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density profile with the kinematic data to the computation of 
the astrophysical factors. This is done using either the GreAT 
MCMC engine or a simpler ;^^/bootstrap approach, depending 
on the user’s choice. 

As for the first release, the code is fully documented and pro¬ 
vides many examples. Many options are available in command 
line executables, and several output formats (ROOT C-r-h based 
plots, fits files and/or ASCII and .root files) are now avail¬ 
able, hopefully making CLUMPY a user-friendly code for both 
the particle physics and astrophysics communities. 
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